home *** CD-ROM | disk | FTP | other *** search
/ Almathera Ten Pack 3: CDPD 3 / Almathera Ten on Ten - Disc 3: CDPD3.iso / fish / 701-725 / 709 / planets / planet.c < prev    next >
C/C++ Source or Header  |  1995-03-18  |  33KB  |  821 lines
  1. /*******************************************************************************
  2. ** PLANET: A program to determine the position of the Sun, Mercury, Venus,    **
  3. **         Mars, Jupiter, Saturn and Uranus.                                  **
  4. **                                                                            **
  5. ** Reference Jean Meeus's book, ``Astronomical Formulae for Calculators''     **
  6. **                                                                            **
  7. **                               Program by                                   **
  8. **                                                                            **
  9. **                          Fred T. Mendenhall                                **
  10. **                           ihnp4!inuxe!fred                                 **
  11. **                                                                            **
  12. ** Copyright 1985 by Fred T. Mendenhall                                       **
  13. ** All rights reserved                                                        **
  14. **                                                                            **
  15. ** (This copyright prevents you from selling the program -                    **
  16. ** the author grants anyone the right to copy and install the program on any  **
  17. ** machine it will run on)                                                    **
  18. ********************************************************************************
  19. ** Modified by Alan Paeth, awpaeth@watcgl, January 1987 for                   **
  20. **   (1) non-interactive mode (uses current GMT)                              **
  21. **   (2) output in simplified Yale format for use with star charter software  **
  22. **   (3) lines 365-380 rewritten as a simplier C expression (April '87)       **
  23. ********************************************************************************
  24. ** Modified by Petri Launiainen, pl@intrin.FI, February 1987 for              **
  25. **   easier LOGFILE definition/modification in Makefile                       **
  26. ********************************************************************************
  27. ** Modified by Alan Paeth, September 1987, for command line flags.            **
  28. **    The program continues to provide reduced Yale output in the old style,  **
  29. **    but new versions of `starchart' accept either.  The old format allows   **
  30. **    representation of magnitudes below -0.99 (useful for planets), but no   **
  31. **    spectral or stellar identification, which is not useful for planets.    **
  32. ********************************************************************************
  33. ** Modified by Jim Cobb, March 1988, to compute moon positions as well.       **
  34. **    Also, to  change slightly the computation of epli, the obliquity        **
  35. **    of the ecliptic.  Formerly epli had a term                              **
  36. **                                                                            **
  37. **                      +0.00256*cos(omeg*rad).                               **
  38. **                                                                            **
  39. **    This term is a compensation for computing the apparent position of the  **
  40. **    sun (Chapter 18 in Meeus).  But for converting ecliptical coordinates   **
  41. **    to right ascension and declination we want the mean value               **
  42. **    (that is, we don't want this term).                                     **
  43. ********************************************************************************
  44. ** Modified by Bob Leivian, Sept 1989, to compute local azimuth and altitude  **
  45. **     given a lat/lon, also fixed time so you can replace a single item      **
  46. **     without having to supply all (i.e., for today at 9:00 put -t 9)        **
  47. **     also ran the code through 'cb' to standardize the indention.           **
  48. **     I split the code in to more managible pieces.                          **
  49. **     Also printed out time in local and UTC as well as JD.                  **
  50. **                                                                            **
  51. ** See "Time Notes" section below if you plan to rewrite the user interface.  **
  52. ********************************************************************************
  53. ** Modified by Andreas Scherer, July 1992, obviously redoing a lot of Bob's   **
  54. **     work, in that the new AMIGA C-compilers are dealing time according to  **
  55. **     the ANSI standard (Jan 1st 1970!).  At this place I included the       **
  56. **     JULDAY routine and some specialities of the SAS/C-compiler.            **
  57. **     Also I cleaned up the source code for better reading.  Doesn't it look **
  58. **     nicer that way?  I didn't use a code beautifier!                       **
  59. **     Most output texts now have their German counterparts.  You have to     **
  60. **     #define EUTSCH, which gives the command line option `-DEUTSCH'.        **
  61. **     I also included some improvements, mostly concerning Cee, but some     **
  62. **     formulas needed a bit of speed up, like most of the polynomials, which **
  63. **     are now evaluated by the Horner algorithm.                             **
  64. *******************************************************************************/
  65.  
  66. #include <stdio.h>
  67. #include <stdlib.h>
  68. #include <math.h>
  69.  
  70. /*******************************************************************************
  71. ** We need the Longitude, the Latitude, and the Time Zone.                    **
  72. ** Default latitude longitude (of Tempe, AZ) -- (-111.94, 33.32, 7 MST).      **
  73. ** Replace with your favorite place.                                          **
  74. ** Now it's Lichtenfels, Bavaria.                                             **
  75. *******************************************************************************/
  76. #define DEFLONG 11.1
  77. #define DEFLAT  50.1
  78. #define DEFZONE (-1.) /* CET */
  79.  
  80. #define RAD 0.01745329251994329651 /* PI/180 */
  81.  
  82. #define JD1900 2415020. /* Julian date for January 1st, 1900, 12h */
  83. #define JD2000 2451545. /* Julian date for January 1st, 2000, 12h */
  84.  
  85. #define JJ100 36525. /* Days in a Julian century */
  86.  
  87. #ifndef PLANETFILE
  88. #  define PLANETFILE "planet.star"
  89. #endif
  90.  
  91. double longitude = DEFLONG;
  92. double latitude  = DEFLAT;
  93. double GHA_Aries;
  94. int    zone;
  95.  
  96. FILE *logfile;
  97. char *progname;
  98.  
  99. /*******************************************************************************
  100. ** Times Notes:                                                               **
  101. **                                                                            **
  102. ** This software would ideally use a generic UNIX system call to which would  **
  103. ** provide day and date info as integers, plus current time west of GMT.      **
  104. ** These values would become the default settings for each of the switches    **
  105. **                                                                            **
  106. **                      -m  -d  -y  -t  -z.                                   **
  107. **                                                                            **
  108. ** Unfortunately, time on across many UNIX systems does not (to my knowledge) **
  109. ** exist as a uniform subroutine call.  (With ANSI they should. Try it!)      **
  110. **                                                                            **
  111. ** Instead, we presume that the very low level "gettimeofday" (GMT in         **
  112. ** seconds, from 1 Jan 1970) exists on all UNIX installations, and use this   **
  113. ** subroutine to use the present time of day when no command line parameters  **
  114. ** are present, as defaults for the current time have not been filled in.     **
  115. **                                                                            **
  116. ** When any flag-style command line parameters are present, hardcoded         **
  117. ** defaults are substitude for -m, -d, -y -t switches, which may then be      **
  118. ** overridden.  The -z flag is filled in by a companion return value from     **
  119. ** "gettimeofday".  This approach makes for poor (on account of hard coding)  **
  120. ** defaults.                                                                  **
  121. **                                                                            **
  122. ** A worthwhile change would be to simplify the command line control          **
  123. ** structure and default mechanism by using a higher level system call, or    **
  124. ** (safer), write the routine to convert Unix GMT into day, month, year, so   **
  125. ** that it lives within this program.  This would allow for more reasonable   **
  126. ** defaults, and still merely one UNIX system call to get GMT time.  Perhaps  **
  127. ** sources from a UNIX system may be studied to do this correctly.            **
  128. **                                                                            **
  129. ** Expertise in "time" on a specific UNIX system is *NOT* what is required    **
  130. ** -- what *IS* required is knowledge and access to enough independent UNIX   **
  131. ** systems to insure (to a high probability) that the code remains portable.  **
  132. ** In other words, only low-level system calls known to work across a large   **
  133. ** range of systems should be employed.                                       **
  134. *******************************************************************************/
  135.  
  136. /*******************************************************************************
  137. ** For getting the current GMT or UTC.                                        **
  138. *******************************************************************************/
  139. #include <time.h>
  140.  
  141. /*******************************************************************************
  142. ** This string is especially for the SAS/C-compiler.  The first three chars   **
  143. ** are the name of your local time zone, the next three represent the normal  **
  144. ** difference from UTC to your local time, and the last three characters      **
  145. ** determine, whether the so called `daylight saving time' is momentarily in  **
  146. ** effect, which means that one hour more is substracted from the UTC value.  **
  147. *******************************************************************************/
  148. char *_TZ = "CET-02CDT";
  149.  
  150. double current_time(m,d,y,t,z)
  151. int    *m,*d,*y,*z;
  152. double *t;
  153.    {
  154.    double julday();
  155.    int    to_int();
  156.  
  157. #ifdef AMIGA
  158. /*******************************************************************************
  159. ** The old version didn't work with the new ANSI-C-compiler. time doesn't     **
  160. ** return the seconds from January 1st, 1978 (happy birthday, AMIGA) but      **
  161. ** the seconds from January 1st, 1970.  We use it here solely for localtime.  **
  162. ** The julday-routine does all the computing.                                 **
  163. ** The first instruction sets up a handful of variables in `time.h', one of   **
  164. ** them is `timezone', the number of seconds, by which UTC is off.  This is   **
  165. ** determined by the `_TZ'-string, initialized above.                         **
  166. *******************************************************************************/
  167.    time_t amiga_now;
  168.    struct tm *p;
  169.  
  170.    tzset();
  171.  
  172.    amiga_now = time(&amiga_now);
  173.    p = gmtime(&amiga_now);
  174.  
  175.    *y = p->tm_year + 1900;          /* local year             */
  176.    *m = p->tm_mon  + 1;             /* local month            */
  177.    *d = p->tm_mday;                 /* local day of month     */
  178.    *t = p->tm_hour + p->tm_min/60.; /* local time (0.-23.999) */
  179.    *z = to_int(timezone / 60.);     /* zone difference in min */
  180.  
  181.    return(julday(*d,*m,*y,*t));
  182.  
  183. #endif /* AMIGA */
  184.  
  185. #ifdef UNIX
  186. /*******************************************************************************
  187. ** I don't know whether the time-routine works in the ANSI-standard.          **
  188. ** It better should. Just in case it doesn't I hold the old version.          **
  189. *******************************************************************************/
  190.    struct timeval tv;
  191.    struct timezone tz;
  192.    struct tm *p;
  193.  
  194.    gettimeofday(&tv,&tz);
  195.    p = localtime(&tv);
  196.  
  197.    *y = p->tm_year + 1900;          /* local year                   */
  198.    *m = p->tm_mon + 1;              /* local month                  */
  199.    *d = p->tm_mday;                 /* local day of month           */
  200.    *t = p->tm_hour + p->tm_min/60.; /* local time (0.-23.999)       */
  201.    *z = tz.tz_minuteswest;          /* local time west of Greenwich */
  202.  
  203.    return(julday(*d,*m,*y,*t + *z));
  204.  
  205. #endif /* UNIX */
  206.  
  207. #ifdef NO_TIME_AVAIL
  208. /*******************************************************************************
  209. ** No time is available, pick a default, i. e., Noon January 1st, 1989        **
  210. ** Due to JULDAY you may use any other set of values, maybe your birthday?    **
  211. *******************************************************************************/
  212.    *y = 1989;                  /* local year                        */
  213.    *m = 1;                     /* local month                       */
  214.    *d = 1;                     /* local day of month                */
  215.    *t = 12.;                   /* local time                        */
  216.    *z = to_int(DEFZONE * 60.); /* zone difference west of Greenwich */
  217.  
  218.    return(julday(*d,*m,*y,*t + DEFZONE);
  219.  
  220. #endif /* NO_TIME_AVAIL */
  221.    }
  222.  
  223. /*******************************************************************************
  224. ** The user is allowed to change the default time settings from his system by **
  225. ** some switches. If some values are changed, they are checked and then they  **
  226. ** may replace the original values.                                           **
  227. *******************************************************************************/
  228. double commandtime(argc,argv)
  229. int argc;
  230. char **argv;
  231.    {
  232.    double time,timez,now;
  233.    int    mm,dd,yy,j;
  234.    double htod(),current_time(),julday();
  235.    int    to_int();
  236.    void   die();
  237.  
  238. #ifdef EUTSCH
  239.    static char *usage =
  240.       "\nAufruf: planet {Zeit}"
  241.       "\noder\tplanet Julianisches Datum"
  242.       "\noder\tplanet [-t hh.mm -z Zone -y Jahr"
  243.       " -m Monat -d Tag -lo Laenge -la Breite]";
  244. #else
  245.    static char *usage =
  246.       "\nusage:planet {current time used}"
  247.       "\nor\tplanet juliandate"
  248.       "\nor\tplanet [-t hh.mm -z zone -y year"
  249.       " -m mon -d day -lo longitude -la latitude]";
  250. #endif
  251.  
  252. /*******************************************************************************
  253. ** Get the current (system) time for defaults.                                **
  254. ** You're getting the fraction of Julian days and the time zone in hours.     **
  255. *******************************************************************************/
  256.    now = current_time(&mm,&dd,&yy,&time,&zone);
  257.    timez = zone / 60.;
  258.  
  259. /*******************************************************************************
  260. ** No args - use current time.                                                **
  261. *******************************************************************************/
  262.    if(argc == 1)
  263.       return(now);
  264.  
  265. /*******************************************************************************
  266. ** One numeric arg - use it as the Julian date.                               **
  267. *******************************************************************************/
  268.    if((argv[1][0] != '-')) {
  269.       if(argc == 2)
  270.          return(atof(argv[1]));
  271.       else
  272. /*******************************************************************************
  273. ** Too many args without switches.                                            **
  274. *******************************************************************************/
  275. #ifdef EUTSCH
  276.          die("Keine Schalter - %s",usage);
  277. #else
  278.          die("no switches - %s",usage);
  279. #endif
  280.       }
  281.  
  282. /*******************************************************************************
  283. ** Flags present, set up defaults, process any user overrides.                **
  284. *******************************************************************************/
  285.    else {
  286.       for(j=1; j<argc; j++) {
  287.          if(argv[j][0] != '-')
  288. #ifdef EUTSCH
  289.             die("Unbekanntes Argument - %s",usage);
  290. #else
  291.             die("unknown argument - %s",usage);
  292. #endif
  293.          switch(argv[j][1]) {
  294. /*******************************************************************************
  295. ** Change the time of day.                                                    **
  296. *******************************************************************************/
  297.          case 't': 
  298.             time = htod(argv[++j]); break;
  299. /*******************************************************************************
  300. ** Change the current time zone.                                              **
  301. *******************************************************************************/
  302.          case 'z': 
  303.             timez = htod(argv[++j]);
  304.             if(timez <= -12 || timez > 12)
  305. #ifdef EUTSCH
  306.                die("Zeitzone: GMT +- 12h (oestlich negativ)\n");
  307. #else
  308.                die("time zone: UTC +- 12h (east negativ)\n");
  309. #endif
  310.             zone = timez * 60.;
  311.             break;
  312. /*******************************************************************************
  313. ** Change the current year.                                                   **
  314. *******************************************************************************/
  315.          case 'y': 
  316.             yy = atoi(argv[++j]); break;
  317. /*******************************************************************************
  318. ** Change the current month of the year.                                      **
  319. *******************************************************************************/
  320.          case 'm': 
  321.             mm = atoi(argv[++j]);
  322.             if(mm < 1 || mm > 12)
  323. #ifdef EUTSCH
  324.                die("Monat zwischen 1 und 12\n");
  325. #else
  326.                die("month must be 1..12\n");
  327. #endif
  328.                break;
  329. /*******************************************************************************
  330. ** Change the current day.                                                    **
  331. *******************************************************************************/
  332.          case 'd': 
  333.             dd = atoi(argv[++j]);
  334.             if(dd < 1 || dd > 31)
  335. #ifdef EUTSCH
  336.                die("Tag zwischen 1 und 31\n");
  337. #else
  338.                die("day must be 1..31\n");
  339. #endif
  340.             break;
  341. /*******************************************************************************
  342. ** Change the position of observation.                                        **
  343. *******************************************************************************/
  344.          case 'l':
  345.             switch(argv[j][2]) {
  346. /*******************************************************************************
  347. ** Change the default longitude of the observer's position.                   **
  348. *******************************************************************************/
  349.             case 'o':
  350.                longitude = atof(argv[++j]);
  351.                if(longitude < -180 || longitude > 180)
  352. #ifdef EUTSCH
  353.                   die("Laenge zwischen -180 und +180\n");
  354. #else
  355.                   die("longitude must be +-180\n");
  356. #endif
  357.                break;
  358. /*******************************************************************************
  359. ** Change the default latitude of the observer's position.                    **
  360. *******************************************************************************/
  361.             case 'a':
  362.                latitude = atof(argv[++j]);
  363.                if(latitude < -90 || latitude > 90)
  364. #ifdef EUTSCH
  365.                   die("Breite zwischen -90 und +90\n");
  366. #else
  367.                   die("latitude must be +-90\n");
  368. #endif
  369.                break;
  370. /*******************************************************************************
  371. ** Tell me what you want.                                                     **
  372. *******************************************************************************/
  373.             default:
  374. #ifdef EUTSCH
  375.                die("Breite und Laenge angeben\n");
  376. #else
  377.                die("specify lat or lon\n");
  378. #endif
  379.                }
  380.             break;
  381. /*******************************************************************************
  382. ** Something went wrong.                                                      **
  383. *******************************************************************************/
  384.          default:
  385. #ifdef EUTSCH
  386.             die("Unbekannter Schalter - %s",usage);
  387. #else
  388.             die("unknown switch - %s",usage);
  389. #endif
  390.             }
  391. /*******************************************************************************
  392. ** ???                                                                        **
  393. *******************************************************************************/
  394.          if(j == argc)
  395. #ifdef EUTSCH
  396.             die("Schalter zuviel - %s",argv[j-1]);
  397. #else
  398.             die("trailing command line flag - %s",argv[j - 1]);
  399. #endif
  400.          }
  401.       }
  402.    return(julday(dd,mm,yy,time + timez));
  403.    }
  404.  
  405. /*******************************************************************************
  406. ** It's much easier to `clear the screen' than most people think.             **
  407. ** There might be some trouble on IBM's, but who cares?                       **
  408. *******************************************************************************/
  409. void cls()
  410.    {
  411.    putchar('\f');
  412.    }
  413.  
  414. /*******************************************************************************
  415. ** Compute the value of a polynomial with degree 3 at time t.                 **
  416. *******************************************************************************/
  417. double poly(a0,a1,a2,a3,t)
  418. double a0,a1,a2,a3,t;
  419.    {
  420.    return(a0 + (a1 + (a2 + a3 * t) * t) * t);
  421.    }
  422.  
  423. /*******************************************************************************
  424. ** Return the integer part of a floating point number.                        **
  425. *******************************************************************************/
  426. int to_int(z)
  427. double z;
  428.    {
  429.    return((int)z);
  430.    }
  431.  
  432. int to_long(z)
  433. double z;
  434.    {
  435.    return((long)z);
  436.    }
  437.  
  438. double kepler(e,M)
  439. double e,M;
  440.    {
  441.    double corr=1.,E0=M,E1;
  442.  
  443.    while(corr > 0.000001) {
  444.       corr = (M + e/RAD * sin(E0*RAD) - E0) / (1 - e * cos(E0*RAD));
  445.       E1 = E0 + corr;
  446.       if(corr < 0)
  447.          corr *= -1.;
  448.       E0 = E1;
  449.       }
  450.    return(E1);
  451.    }
  452.  
  453. double truean(e,E)
  454. double e,E;
  455.    {
  456.    double nu;
  457.  
  458.    nu  = 2. * atan(sqrt((1 + e) / (1 - e)) * tan(E * RAD / 2)) / RAD;
  459.    if(nu < 0.)
  460.       nu += 360.;
  461.    if(fabs(nu - E) > 90.)
  462.       if(nu > 180.)
  463.          nu -= 180.;
  464.       else
  465.          nu += 180.;
  466.    return(nu);
  467.    }
  468.  
  469. double longi(w2,i,u)
  470. double w2,i,u;
  471.    {
  472.    double x,y,l;
  473.  
  474.    y = cos(i*RAD) * sin(u*RAD);
  475.    x = cos(u*RAD);
  476.    l = atan2(y,x)/RAD;
  477.    if(l<0.)
  478.       l += 360.;
  479.    return(l + w2);
  480.    }
  481.  
  482. double lati(u,i)
  483. double u,i;
  484.    {
  485.    return(asin(sin(u*RAD) * sin(i*RAD)) / RAD);
  486.    }
  487.  
  488. void speak(ra,dec,dis,mag,sym,str)
  489. double ra,dec,dis,mag;
  490. char *sym,*str;
  491.    {
  492.    double lha,sha,altitude,azimuth;
  493.    int    h,m,s,deg;
  494.    double range();
  495.    int    to_int();
  496.  
  497.    if(ra < 0)
  498.       ra += 360.;
  499.    sha = 360. - ra;
  500.  
  501. /*******************************************************************************
  502. ** Convert from degs to hours.                                                **
  503. *******************************************************************************/
  504.    ra /= 15.;
  505.    h = to_int(ra);
  506.    m = to_int( (ra - h) * 60          );
  507.    s = to_int(((ra - h) * 60 - m) * 60);
  508.  
  509.    printf(" %-12s%2dh %02dm %02ds ",str,h,m,s);
  510.  
  511.    if(logfile)
  512.       fprintf(logfile,"%02d%02d%02d",h,m,s);
  513.  
  514.    deg = to_int(dec);
  515.    m   = to_int( (dec - deg) * 60          );
  516.    s   = to_int(((dec - deg) * 60 - m) * 60);
  517.    if(m < 0)
  518.       m *= -1;
  519.    if(s < 0)
  520.       s *= -1;
  521.  
  522. #ifdef AMIGA
  523. /*******************************************************************************
  524. ** Its printf doesn't know about signed. Well, I didn't try it yet.           **
  525. ***************************************************************************** */
  526.    if(deg > 0)
  527.       printf("   +%2ddeg %02dm %02ds ",deg,m,s);
  528.    else
  529.       printf("   %3ddeg %02dm %02ds ",deg,m,s);
  530. #else
  531.    printf("   %+3ddeg %02dm %02ds ",deg,m,s);
  532. #endif /* AMIGA */
  533.  
  534.    if(dis > 0.0001)
  535.       printf("  %6.3f",dis);
  536.    else
  537.       printf("        ");
  538.    if(logfile)
  539.       if(mag < 0.)
  540.          fprintf(logfile,"%+03d%02d-%2d%s       %s\n",
  541.             deg,m,-to_int(10 * mag),sym,str);
  542.       else
  543.          fprintf(logfile,"%+03d%02d%3d%s       %s\n",
  544.             deg,m,to_int(100 * mag),sym,str);
  545.  
  546. /*******************************************************************************
  547. ** If we have a valid angle.                                                  **
  548. *******************************************************************************/
  549.    if(GHA_Aries > 0) {
  550.  
  551. /*******************************************************************************
  552. ** Now tell where it is at a given place on earth.                            **
  553. *******************************************************************************/
  554.       lha = range(GHA_Aries + sha + longitude);
  555.  
  556.       altitude = sin(latitude * RAD) * sin(dec * RAD) + 
  557.          cos(latitude * RAD) * cos(dec * RAD) * cos(lha * RAD);
  558.       altitude = asin(altitude) / RAD;
  559.  
  560.       azimuth = sin(lha * RAD) / 
  561.          (cos(lha * RAD) * sin(latitude * RAD) - 
  562.          tan(dec * RAD) * cos(latitude * RAD));
  563.  
  564. /*******************************************************************************
  565. ** Correct for quadrant.                                                      **
  566. *******************************************************************************/
  567.       if(lha <= 180.) {
  568.          if(azimuth > 0)
  569.             azimuth = 180. + atan(azimuth) / RAD;
  570.          else
  571.             azimuth = 360. + atan(azimuth) / RAD;
  572.          }
  573.       else {
  574.          if(azimuth <= 0.)
  575.             azimuth = 180. + atan(azimuth) / RAD;
  576.          else
  577.             azimuth = atan(azimuth) / RAD;
  578.          }
  579.  
  580.       printf("      %5.1f     %5.1f\n",altitude,azimuth);
  581.       }
  582.    else
  583.       printf("\n");
  584.    }
  585.  
  586. /*******************************************************************************
  587. ** Helio_trans converts heliocentric ecliptical coordinates to geocentric     **
  588. ** right ascension and declination.  It also outputs the converted value.     **
  589. *******************************************************************************/
  590. void helio_trans(r,b,ll,Stheta,Sr,epli,mag,sym,str)
  591. double r,b,ll,Stheta,Sr,epli,mag;
  592. char   *sym,*str;
  593.    {
  594.    double N,D,lambda,delta,beta,rcosb=r*cos(b*RAD),rsinb=r*sin(b*RAD);
  595.    double range();
  596.    void   geo_trans();
  597.  
  598.    N = rcosb * sin((ll - Stheta)*RAD);
  599.    D = rcosb * cos((ll - Stheta)*RAD) + Sr;
  600.    lambda = atan2(N,D) / RAD;
  601.    if(lambda<0.)
  602.       lambda += 360.;
  603.    lambda = range(lambda + Stheta);
  604.    delta  = sqrt(N*N + D*D + rsinb * rsinb);
  605.    beta   = asin(rsinb / delta) / RAD;
  606.    geo_trans(lambda,beta,epli,delta,mag,sym,str);
  607.    }
  608.  
  609. /*******************************************************************************
  610. ** geo_trans converts geocentric ecliptical coordinates to geocentric right   **
  611. ** ascension and declination.  It also outputs the converted value.  Note     **
  612. ** that the output coordinates are referred to the mean equinox of date.  If  **
  613. ** these coordinates are to be used in conjunction with a star database, use  **
  614. ** the epoch program to convert the planet's coordinates to the epoch for the **
  615. ** star database.                                                             **
  616. *******************************************************************************/
  617. void geo_trans(lambda,beta,epli,delta,mag,sym,str)
  618. double lambda,beta,epli,delta,mag;
  619. char   *sym,*str;
  620.    {
  621.    double N,D,RA,DEC;
  622.    double sinlam=sin(lambda*RAD),cosepl=cos(epli*RAD),sinepl=sin(epli*RAD);
  623.    void speak();
  624.  
  625.    N   = sinlam * cosepl - tan(beta*RAD) * sinepl;
  626.    D   = cos(lambda*RAD);
  627.    RA  = atan2(N,D)/RAD;
  628.    DEC = asin(sin(beta*RAD)*cosepl + cos(beta*RAD)*sinepl*sinlam)/RAD;
  629.    speak(RA,DEC,delta,mag,sym,str);
  630.    }
  631.  
  632. /*******************************************************************************
  633. ** htod(str) reads floating point strings of the form {+-}hh.{mm} thus        **
  634. ** allowing for fractional values in base 60 (i.e., degrees/minutes).         **
  635. *******************************************************************************/
  636. double htod(s)
  637. char *s;
  638.    {
  639.    double sign;
  640.    int    full,frac;
  641.    char   *t;
  642.    void   die();
  643.  
  644.    t = s-1;
  645.    while(*++t) {
  646.       if(((*t<'0') || (*t>'9')) && (*t!='.') && (*t!='+') && (*t!='-'))
  647. #ifdef EUTSCH
  648.          die("Falsches Zeichen im numerischen Argument tt.mm");
  649. #else
  650.          die("non-digit in dd.mm style numeric argument");
  651. #endif
  652.       }
  653.    if(s == NULL)
  654.       return(0.);
  655.  
  656.    full = frac = 0;
  657.    sign = 1.;
  658.  
  659.    if(*s == '-') {
  660.       sign = -1.;
  661.       s++;
  662.       }
  663.    else if(*s == '+')
  664.       s++;
  665.  
  666.    while(*s && *s != '.')
  667.       full = 10*full + *s++ - '0';
  668.  
  669.    if(*s++ == '.') {
  670.       if(*s)
  671.          frac = 10*(*s++ - '0');
  672.       if(*s)
  673.          frac += *s++ - '0';
  674.       if(frac>59)
  675.          frac = 59;
  676.       }
  677.    return(sign * ((double)full + ((double)frac) / 60.));
  678.    }
  679.  
  680. /*******************************************************************************
  681. ** In case of emergency you need a safety exit.                               **
  682. *******************************************************************************/
  683. void die(a,b)
  684. char *a,*b;
  685.    {
  686.    fprintf(stderr,"%s: ",progname);
  687.    fprintf(stderr,a,b);
  688.    fprintf(stderr,"\n");
  689.    exit(1);
  690.    }
  691.  
  692. /*******************************************************************************
  693. ** The main routine.  I placed it to the end.                                 **
  694. *******************************************************************************/
  695. void main(argc,argv)
  696. int argc;
  697. char **argv;
  698.    {
  699.    double jd,T0,epli,temp;
  700.    long   int_jd;
  701.    int    y,m,d,h,min,lh,lmin;
  702.    char   *str;
  703.    double planet_pos(),sidereal_time(),commandtime();
  704.    int    to_int();
  705.    void   cls(),speak(),die(),moon_pos(),caldat(),caltim();
  706.  
  707.    progname = argv[0];
  708.  
  709.    if((logfile = fopen(PLANETFILE,"w")) == 0)
  710. #ifdef EUTSCH
  711.       die("Fehler beim Oeffnen der Logdatei %s\n",PLANETFILE);
  712. #else
  713.       die("fail on opening logfile %s\n",PLANETFILE);
  714. #endif
  715.  
  716. /*******************************************************************************
  717. ** Calculate Julian date in fractions of days.                                **
  718. *******************************************************************************/
  719.    jd = commandtime(argc,argv);
  720.  
  721. /*******************************************************************************
  722. ** Convert the Julian date to calendar day.                                   **
  723. *******************************************************************************/
  724.    int_jd = to_int(jd +.5001);
  725.    caldat(int_jd,&m,&d,&y);
  726.  
  727. /*******************************************************************************
  728. ** Julian day to UTC.                                                         **
  729. *******************************************************************************/
  730.    temp = jd - .4999;
  731.    caltim(temp - floor(temp),&h,&min);
  732.  
  733.    cls();
  734.  
  735. #ifdef EUTSCH
  736.    printf("\n Am %d.%02d.%4d um %d Uhr %02d GMT (Julianisch: %.3f)\n",
  737.       d,m,y,h,min,jd);
  738. #else
  739.    printf("\n At %d/%d %d %d:%02d UTC (Julian: %.3f)\n",
  740.       m,d,y,h,min,jd);
  741. #endif
  742.  
  743. /*******************************************************************************
  744. ** UTC to Local.                                                              **
  745. *******************************************************************************/
  746.    temp -= (zone / 1440.);
  747.    caltim(temp - floor(temp),&lh,&lmin);
  748.    int_jd = to_int(jd + .5001 - (zone / 1440.));
  749.    caldat(int_jd,&m,&d,&y);
  750.  
  751. #ifdef EUTSCH
  752.    str = "keine Zeitangabe";
  753. #else
  754.    if(lh > 12) {
  755.       lh -= 12;
  756.       str = "pm";
  757.       }
  758.    else {
  759.       str = "am";
  760.       if(lh == 0)
  761.          lh = 12;
  762.       }
  763. #endif
  764.  
  765. #ifdef EUTSCH
  766.    printf(" Lokal %d.%02d. %d Uhr %02d (",
  767.       d,m,lh,lmin);
  768. #else
  769.    printf(" local %d/%d %d:%02d%s (",
  770.       m,d,lh,lmin,str);
  771. #endif
  772.  
  773. /*******************************************************************************
  774. ** Compute exact local time for a given longitude.                            **
  775. *******************************************************************************/
  776.    temp = jd - 0.5 + (longitude / 360.);
  777.    caltim(temp - floor(temp),&lh,&lmin);
  778.    GHA_Aries = 15. * sidereal_time(jd);
  779.  
  780. #ifdef EUTSCH
  781.    printf("%d:%02d wahr lokal bei Brt: %1.1f Lng: %1.1f, GHA Aries: %1.1f)\n",
  782.       lh,lmin,latitude,longitude,GHA_Aries);
  783. #else
  784.    printf("%d:%02d true local at lat:%1.1f lon:%1.1f, GHA Aries:%1.1f)\n",
  785.       lh,lmin,latitude,longitude,GHA_Aries);
  786. #endif
  787.  
  788.    T0 = (jd - JD1900) / JJ100;
  789.  
  790. #ifdef EUTSCH
  791.    printf("\n Objekt     Rektaszension     Deklination    Abstand     Höhe     Azimuth\n");
  792. #else
  793.    printf("\n OBJECT          RA           DEC           DISTANCE   altitude   azimuth\n");
  794. #endif
  795.  
  796.    epli = planet_pos(jd,T0);
  797.  
  798.    moon_pos(T0,epli);
  799.  
  800.    putchar('\n');
  801.    if(logfile)
  802.       close(logfile);
  803.    logfile = 0;
  804.  
  805. /*******************************************************************************
  806. ** Show local position of some stars also.                                    **
  807. *******************************************************************************/
  808.    speak(101.1714,-16.7013,0.,-1.46," ","Sirius");
  809. #ifdef EUTSCH
  810.    speak(213.7956,19.2360,0.,-0.04," ","Arktur");
  811.    speak( 88.6508, 7.4057,0., 0.50," ","Beteigeuze");
  812. #else
  813.    speak(213.7956,19.2360,0.,-0.04," ","Arcturus");
  814.    speak( 88.6508, 7.4057,0., 0.50," ","Betelgeuse");
  815. #endif
  816.    putchar('\n');
  817.    }
  818. /*******************************************************************************
  819. ** End of the main program and end of this file.                              **
  820. *******************************************************************************/
  821.